library("tidyr")
library("ICSNP")
Loading required package: mvtnorm
Loading required package: ICS
library("robustbase")
library("ggplot2")
SIMULATED DATA
library("pracma")
library("MASS")
# MODEL 1
sepa = 50
t = linspace(0,1,sepa)
dift = max(diff(t))
c = 0.2
n = 100
clean = ceiling(n*(1-c))
outs = n - clean
noise1 = 0.3*exp(-abs(dift)/0.3)
mu = rep(0,sepa)
sigma = matrix(noise1, sepa, sepa)
for (i in 1:sepa) {
for (j in 1:sepa) {
sigma[i,j] = 0.3*exp(-abs(t[i]-t[j])/0.3)
}
}
Main1 = matrix(0, nrow = n, ncol = length(t))
for (i in 1:clean) {
Main1[i,] = 30*t*((1-t)^(3/2)) + mvrnorm(1, mu, sigma)
}
for (i in (clean+1):n) {
Main1[i,] = 30*(t^(3/2))*(1-t) + mvrnorm(1, mu, sigma)
}
#plot(Main1[1,], type="l", col="black", lwd=2)
#lines(-Main1[2,], col="red", lwd=2)
#lines(Main1[3,], col="blue", lwd=2)
data <- c()
for (i in (clean+1):n) {
if (mod(i, 2) == 0) t = 1
else t=-1
data <- c(data, t*Main1[i,])
}
plot(data, type="l", lwd=2)

# MODEL 2
sepa = 50
t = linspace(0,1,sepa)
dift = max(diff(t))
c = 0.2
n = 100
clean = ceiling(n*(1-c))
outs = n - clean
noise1 = exp(-abs(dift))
mu = rep(0,sepa)
sigma = matrix(noise1, sepa, sepa)
for (i in 1:sepa) {
for (j in 1:sepa) {
sigma[i,j] = exp(-abs(t[i]-t[j]))
}
}
Main1 = matrix(0, nrow = n, ncol = length(t))
for (i in 1:clean) {
Main1[i,] = 4*t + mvrnorm(1, mu, sigma)
}
for (i in (clean+1):n) {
Main1[i,] = 4*t + 1.8*((-1)^rbinom(1,1,0.5)) + (1/sqrt(2*pi*0.01))*exp(-((t-runif(1,0.25,0.75))^2)/0.02) + mvrnorm(1, mu, sigma)
}
data <- c()
for (i in (clean+1):n) {
if (mod(i, 2) == 0) t = 1
else t=-1
data <- c(data, t*Main1[i,])
}
plot(data, type="l", lwd=2)

# MODEL 3
sepa = 50
t = linspace(0,1,sepa)
dift = max(diff(t))
c = 0.2
n = 100
clean = ceiling(n*(1-c))
outs = n - clean
noise1 = exp(-abs(dift))
mu = rep(0,sepa)
sigma = matrix(noise1, sepa, sepa)
for (i in 1:sepa) {
for (j in 1:sepa) {
sigma[i,j] = exp(-abs(t[i]-t[j]))
}
}
Main1 = matrix(0, nrow = n, ncol = length(t))
for (i in 1:clean) {
Main1[i,] = 4*t + mvrnorm(1, mu, sigma)
}
for (i in (clean+1):n) {
Main1[i,] = 4*t + 2*sin(4*(t+runif(1,0.25,0.75))*pi) + mvrnorm(1, mu, sigma)
}
data <- c()
for (i in (clean+1):n) {
if (mod(i, 2) == 0) t = 1
else t=-1
data <- c(data, t*Main1[i,])
}
plot(data, type="l", lwd=2)

data_cont <- data
scales <- matrix(ncol=6)
for (i in seq(0.01, 0.4, 0.01)) {
alfa <- round(length(data) * i)
for (j in seq(1, alfa, 1)){
w = round(runif(1, 1, length(data)))
if (mod(w, 2) == 0) data_cont[w] <- data_cont[w]*1.15
#else data_cont[w]<- data_cont[w]*-sd(data)/10
}
scales <- rbind(scales, c(Pn(data_cont), Qn(data_cont), Sn(data_cont), mad(data_cont), sd(data_cont), IQR(data_cont)))
}
scales <- cbind(scales, seq(1, nrow(scales)))
colnames(scales) <- c("Pn", "Qn", "Sn", "MAD", "SD", "IQR", "X")
df <- data.frame(scales) %>%
gather(key, value, -X)
p <- ggplot(df, aes(X, value, color = key)) + geom_line()
plotly::ggplotly(p)
#plot(data_cont, type="l", lwd=1)
plot(data_cont, type="l", lwd=1)

ROBUST SCALE ESTIMATORS COMPARISON Pn extracted from https://github.com/garthtarr/robnetwork/blob/master/covest.R using Pn scale estimator
Pn = function(y){
n = length(y)
y.pairs = pair.sum(matrix(y))/2
c_t = 1/0.9539
scale_estim = c_t*as.numeric(diff(quantile(y.pairs,c(0.25,0.75),type=1)))
correction.factors =
c(1.128,1.303,1.109,1.064,1.166,1.103,1.087,1.105,1.047,1.063,1.057,1.040,
1.061,1.047,1.043,1.048,1.031,1.037,1.035,1.028,1.036,1.030,1.029,1.032,
1.023,1.025,1.024,1.021,1.026,1.022,1.021,1.023,1.018,1.020,1.019,1.017,
1.020,1.018,1.017,1.018,1.015,1.016,1.016,1.014,1.016,1.015,1.014,1.015)
if(n <= 40){
scale_estim = scale_estim*correction.factors[n-2]
} else if(n > 40) scale_estim = scale_estim*n/(n-0.7)
return(scale_estim)
}
print(Pn(data))
print(Qn(data))
print(Sn(data))
print(mad(data))
N <- 10000
mu <- 2
sigma <- 1/2
scales <- matrix(ncol=6)
pn <- c()
X = rnorm(N, mean=mu, sd = sigma)
for (i in seq(0.01, 0.2, by=0.01)){
Y = rexp(N*i, rate=mu)
X_Y = c(X, Y)
pn <- c(pn, Pn(X_Y))
scales <- rbind(scales, c(Pn(X_Y), Qn(X_Y), Sn(X_Y), mad(X_Y), sd(X_Y), IQR(X_Y)))
}
scales <- cbind(scales, seq(1, 21))
colnames(scales) <- c("Pn", "Qn", "Sn", "MAD", "SD", "IQR", "X")
df <- data.frame(scales) %>%
gather(key, value, -X)
p <- ggplot(df, aes(X, value, color = key)) + geom_line()
plotly::ggplotly(p)
N <- 10000
mu <- 2
sigma <- 1/2
scales <- matrix(ncol=6)
pn <- c()
#X = rnorm(N, mean=mu, sd = sigma)
X_w = rweibull(N, shape=mu, scale=sigma)
for (i in seq(0.01, 0.2, by=0.01)){
Y = rexp(N*i, rate=mu)
X_Y = c(X_w, Y)
pn <- c(pn, Pn(X_Y))
scales <- rbind(scales, c(Pn(X_Y), Qn(X_Y), Sn(X_Y), mad(X_Y), sd(X_Y), IQR(X_Y)))
}
scales <- cbind(scales, seq(1, 21))
colnames(scales) <- c("Pn", "Qn", "Sn", "MAD", "SD", "IQR", "X")
df <- data.frame(scales) %>%
gather(key, value, -X)
p <- ggplot(df, aes(X, value, color = key)) + geom_line()
plotly::ggplotly(p)
N <- 10000
mu <- 2
sigma <- 1/2
scales <- matrix(ncol=6)
pn <- c()
#X = rnorm(N, mean=mu, sd = sigma)
X_w = rweibull(N, shape=mu, scale=sigma)
for (i in seq(0.01, 0.2, by=0.01)){
Y = rcauchy(N*i, location=mu, scale=sigma/1000)
X_Y = c(X_w, Y)
pn <- c(pn, Pn(X_Y))
scales <- rbind(scales, c(Pn(X_Y), Qn(X_Y), Sn(X_Y), mad(X_Y), sd(X_Y), IQR(X_Y)))
}
scales <- cbind(scales, seq(1, 21))
colnames(scales) <- c("Pn", "Qn", "Sn", "MAD", "SD", "IQR", "X")
df <- data.frame(scales) %>%
gather(key, value, -X)
p <- ggplot(df, aes(X, value, color = key)) + geom_line()
plotly::ggplotly(p)
LS0tDQp0aXRsZTogIlNjYWxlcyBFc3RpbWF0b3JzIHZzIENvbnRhbWluYXRpb24iDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpgYGB7cn0NCmxpYnJhcnkoInRpZHlyIikNCmxpYnJhcnkoIklDU05QIikNCmxpYnJhcnkoInJvYnVzdGJhc2UiKQ0KbGlicmFyeSgiZ2dwbG90MiIpDQpsaWJyYXJ5KCJwcmFjbWEiKQ0KbGlicmFyeSgiTUFTUyIpDQpgYGANCg0KU0lNVUxBVEVEIERBVEENCg0KDQpgYGB7cn0NCiMgTU9ERUwgMQ0Kc2VwYSA9IDUwDQp0ID0gbGluc3BhY2UoMCwxLHNlcGEpDQpkaWZ0ID0gbWF4KGRpZmYodCkpDQpjID0gMC4yDQpuID0gMTAwDQpjbGVhbiA9IGNlaWxpbmcobiooMS1jKSkNCm91dHMgPSAgbiAtIGNsZWFuDQpub2lzZTEgPSAwLjMqZXhwKC1hYnMoZGlmdCkvMC4zKQ0KbXUgPSByZXAoMCxzZXBhKQ0Kc2lnbWEgPSBtYXRyaXgobm9pc2UxLCBzZXBhLCBzZXBhKQ0KZm9yIChpIGluIDE6c2VwYSkgew0KICBmb3IgKGogaW4gMTpzZXBhKSB7DQogICAgc2lnbWFbaSxqXSA9IDAuMypleHAoLWFicyh0W2ldLXRbal0pLzAuMykNCiAgfQ0KfQ0KTWFpbjEgPSBtYXRyaXgoMCwgbnJvdyA9IG4sIG5jb2wgPSBsZW5ndGgodCkpDQpmb3IgKGkgaW4gMTpjbGVhbikgew0KICBNYWluMVtpLF0gPSAzMCp0KigoMS10KV4oMy8yKSkgKyBtdnJub3JtKDEsIG11LCBzaWdtYSkNCn0NCmZvciAoaSBpbiAoY2xlYW4rMSk6bikgew0KICBNYWluMVtpLF0gPSAzMCoodF4oMy8yKSkqKDEtdCkgKyBtdnJub3JtKDEsIG11LCBzaWdtYSkNCn0NCg0KI3Bsb3QoTWFpbjFbMSxdLCB0eXBlPSJsIiwgY29sPSJibGFjayIsIGx3ZD0yKQ0KI2xpbmVzKC1NYWluMVsyLF0sIGNvbD0icmVkIiwgbHdkPTIpDQojbGluZXMoTWFpbjFbMyxdLCBjb2w9ImJsdWUiLCBsd2Q9MikNCmRhdGEgPC0gYygpDQpmb3IgKGkgaW4gKGNsZWFuKzEpOm4pIHsNCiAgaWYgKG1vZChpLCAyKSA9PSAwKSB0ID0gMQ0KICBlbHNlIHQ9LTENCiAgZGF0YSA8LSBjKGRhdGEsIHQqTWFpbjFbaSxdKQ0KfSANCnBsb3QoZGF0YSwgdHlwZT0ibCIsIGx3ZD0yKQ0KYGBgDQoNCg0KDQpgYGB7cn0NCiMgTU9ERUwgMg0Kc2VwYSA9IDUwDQp0ID0gbGluc3BhY2UoMCwxLHNlcGEpDQpkaWZ0ID0gbWF4KGRpZmYodCkpDQpjID0gMC4yDQpuID0gMTAwDQpjbGVhbiA9IGNlaWxpbmcobiooMS1jKSkNCm91dHMgPSAgbiAtIGNsZWFuDQpub2lzZTEgPSBleHAoLWFicyhkaWZ0KSkNCm11ID0gcmVwKDAsc2VwYSkNCnNpZ21hID0gbWF0cml4KG5vaXNlMSwgc2VwYSwgc2VwYSkNCmZvciAoaSBpbiAxOnNlcGEpIHsNCiAgZm9yIChqIGluIDE6c2VwYSkgew0KICAgIHNpZ21hW2ksal0gPSBleHAoLWFicyh0W2ldLXRbal0pKQ0KICB9DQp9DQpNYWluMSA9IG1hdHJpeCgwLCBucm93ID0gbiwgbmNvbCA9IGxlbmd0aCh0KSkNCmZvciAoaSBpbiAxOmNsZWFuKSB7DQogIE1haW4xW2ksXSA9IDQqdCArIG12cm5vcm0oMSwgbXUsIHNpZ21hKQ0KfQ0KZm9yIChpIGluIChjbGVhbisxKTpuKSB7DQogIE1haW4xW2ksXSA9IDQqdCArIDEuOCooKC0xKV5yYmlub20oMSwxLDAuNSkpICsgKDEvc3FydCgyKnBpKjAuMDEpKSpleHAoLSgodC1ydW5pZigxLDAuMjUsMC43NSkpXjIpLzAuMDIpICsgbXZybm9ybSgxLCBtdSwgc2lnbWEpDQp9DQoNCmRhdGEgPC0gYygpDQpmb3IgKGkgaW4gKGNsZWFuKzEpOm4pIHsNCiAgaWYgKG1vZChpLCAyKSA9PSAwKSB0ID0gMQ0KICBlbHNlIHQ9LTENCiAgZGF0YSA8LSBjKGRhdGEsIHQqTWFpbjFbaSxdKQ0KfSANCnBsb3QoZGF0YSwgdHlwZT0ibCIsIGx3ZD0yKQ0KYGBgDQoNCmBgYHtyfQ0KIyBNT0RFTCAzDQpzZXBhID0gNTANCnQgPSBsaW5zcGFjZSgwLDEsc2VwYSkNCmRpZnQgPSBtYXgoZGlmZih0KSkNCmMgPSAwLjINCm4gPSAxMDANCmNsZWFuID0gY2VpbGluZyhuKigxLWMpKQ0Kb3V0cyA9ICBuIC0gY2xlYW4NCm5vaXNlMSA9IGV4cCgtYWJzKGRpZnQpKQ0KbXUgPSByZXAoMCxzZXBhKQ0Kc2lnbWEgPSBtYXRyaXgobm9pc2UxLCBzZXBhLCBzZXBhKQ0KZm9yIChpIGluIDE6c2VwYSkgew0KICBmb3IgKGogaW4gMTpzZXBhKSB7DQogICAgc2lnbWFbaSxqXSA9IGV4cCgtYWJzKHRbaV0tdFtqXSkpDQogIH0NCn0NCk1haW4xID0gbWF0cml4KDAsIG5yb3cgPSBuLCBuY29sID0gbGVuZ3RoKHQpKQ0KZm9yIChpIGluIDE6Y2xlYW4pIHsNCiAgTWFpbjFbaSxdID0gNCp0ICsgbXZybm9ybSgxLCBtdSwgc2lnbWEpDQp9DQpmb3IgKGkgaW4gKGNsZWFuKzEpOm4pIHsNCiAgTWFpbjFbaSxdID0gNCp0ICsgMipzaW4oNCoodCtydW5pZigxLDAuMjUsMC43NSkpKnBpKSArIG12cm5vcm0oMSwgbXUsIHNpZ21hKQ0KfQ0KZGF0YSA8LSBjKCkNCmZvciAoaSBpbiAoY2xlYW4rMSk6bikgew0KICBpZiAobW9kKGksIDIpID09IDApIHQgPSAxDQogIGVsc2UgdD0tMQ0KICBkYXRhIDwtIGMoZGF0YSwgdCpNYWluMVtpLF0pDQp9IA0KcGxvdChkYXRhLCB0eXBlPSJsIiwgbHdkPTIpDQpgYGANCg0KYGBge3J9DQpkYXRhX2NvbnQgPC0gZGF0YQ0Kc2NhbGVzIDwtIG1hdHJpeChuY29sPTYpDQpmb3IgKGkgaW4gc2VxKDAuMDEsIDAuNCwgMC4wMSkpIHsNCiAgYWxmYSA8LSByb3VuZChsZW5ndGgoZGF0YSkgKiBpKQ0KICBmb3IgKGogaW4gc2VxKDEsIGFsZmEsIDEpKXsNCiAgICB3ID0gcm91bmQocnVuaWYoMSwgMSwgbGVuZ3RoKGRhdGEpKSkNCiAgICBpZiAobW9kKHcsIDIpID09IDApIGRhdGFfY29udFt3XSA8LSBkYXRhX2NvbnRbd10qMS4xNQ0KICAgICNlbHNlIGRhdGFfY29udFt3XTwtIGRhdGFfY29udFt3XSotc2QoZGF0YSkvMTANCiAgfQ0KICBzY2FsZXMgPC0gcmJpbmQoc2NhbGVzLCBjKFBuKGRhdGFfY29udCksIFFuKGRhdGFfY29udCksIFNuKGRhdGFfY29udCksIG1hZChkYXRhX2NvbnQpLCBzZChkYXRhX2NvbnQpLCBJUVIoZGF0YV9jb250KSkpDQp9DQpzY2FsZXMgPC0gY2JpbmQoc2NhbGVzLCBzZXEoMSwgbnJvdyhzY2FsZXMpKSkNCmNvbG5hbWVzKHNjYWxlcykgPC0gYygiUG4iLCAiUW4iLCAiU24iLCAiTUFEIiwgIlNEIiwgIklRUiIsICJYIikNCg0KZGYgPC0gZGF0YS5mcmFtZShzY2FsZXMpICU+JQ0KICBnYXRoZXIoa2V5LCB2YWx1ZSwgLVgpDQpwIDwtIGdncGxvdChkZiwgYWVzKFgsIHZhbHVlLCBjb2xvciA9IGtleSkpICsgZ2VvbV9saW5lKCkNCnBsb3RseTo6Z2dwbG90bHkocCkNCiNwbG90KGRhdGFfY29udCwgdHlwZT0ibCIsIGx3ZD0xKQ0KYGBgDQoNCmBgYHtyfQ0KcGxvdChkYXRhX2NvbnQsIHR5cGU9ImwiLCBsd2Q9MSkNCmBgYA0KDQpST0JVU1QgU0NBTEUgRVNUSU1BVE9SUyBDT01QQVJJU09OIA0KUG4gZXh0cmFjdGVkIGZyb20gaHR0cHM6Ly9naXRodWIuY29tL2dhcnRodGFyci9yb2JuZXR3b3JrL2Jsb2IvbWFzdGVyL2NvdmVzdC5SIHVzaW5nIFBuIHNjYWxlIGVzdGltYXRvcg0KDQpgYGB7cn0NClBuID0gZnVuY3Rpb24oeSl7DQogIG4gPSBsZW5ndGgoeSkNCiAgeS5wYWlycyA9IHBhaXIuc3VtKG1hdHJpeCh5KSkvMg0KICBjX3QgPSAxLzAuOTUzOQ0KICBzY2FsZV9lc3RpbSA9IGNfdCphcy5udW1lcmljKGRpZmYocXVhbnRpbGUoeS5wYWlycyxjKDAuMjUsMC43NSksdHlwZT0xKSkpDQogIGNvcnJlY3Rpb24uZmFjdG9ycyA9DQogICAgYygxLjEyOCwxLjMwMywxLjEwOSwxLjA2NCwxLjE2NiwxLjEwMywxLjA4NywxLjEwNSwxLjA0NywxLjA2MywxLjA1NywxLjA0MCwgDQogICAgICAxLjA2MSwxLjA0NywxLjA0MywxLjA0OCwxLjAzMSwxLjAzNywxLjAzNSwxLjAyOCwxLjAzNiwxLjAzMCwxLjAyOSwxLjAzMiwgDQogICAgICAxLjAyMywxLjAyNSwxLjAyNCwxLjAyMSwxLjAyNiwxLjAyMiwxLjAyMSwxLjAyMywxLjAxOCwxLjAyMCwxLjAxOSwxLjAxNywgDQogICAgICAxLjAyMCwxLjAxOCwxLjAxNywxLjAxOCwxLjAxNSwxLjAxNiwxLjAxNiwxLjAxNCwxLjAxNiwxLjAxNSwxLjAxNCwxLjAxNSkNCiAgaWYobiA8PSA0MCl7DQogICAgc2NhbGVfZXN0aW0gPSBzY2FsZV9lc3RpbSpjb3JyZWN0aW9uLmZhY3RvcnNbbi0yXQ0KICB9IGVsc2UgaWYobiA+IDQwKSBzY2FsZV9lc3RpbSA9IHNjYWxlX2VzdGltKm4vKG4tMC43KQ0KICByZXR1cm4oc2NhbGVfZXN0aW0pDQp9DQpgYGANCg0KYGBge3J9DQpwcmludChQbihkYXRhKSkNCnByaW50KFFuKGRhdGEpKQ0KcHJpbnQoU24oZGF0YSkpDQpwcmludChtYWQoZGF0YSkpDQpgYGANCg0KDQpgYGB7cn0NCk4gPC0gMTAwMDANCm11IDwtIDINCnNpZ21hIDwtIDEvMg0Kc2NhbGVzIDwtIG1hdHJpeChuY29sPTYpDQpwbiA8LSBjKCkNClggPSBybm9ybShOLCBtZWFuPW11LCBzZCA9IHNpZ21hKQ0KZm9yIChpIGluIHNlcSgwLjAxLCAwLjIsIGJ5PTAuMDEpKXsNCiAgWSA9IHJleHAoTippLCByYXRlPW11KQ0KICBYX1kgPSBjKFgsIFkpDQogIHBuIDwtIGMocG4sIFBuKFhfWSkpDQogIHNjYWxlcyA8LSByYmluZChzY2FsZXMsIGMoUG4oWF9ZKSwgUW4oWF9ZKSwgU24oWF9ZKSwgbWFkKFhfWSksIHNkKFhfWSksIElRUihYX1kpKSkNCn0NCnNjYWxlcyA8LSBjYmluZChzY2FsZXMsIHNlcSgxLCAyMSkpDQpjb2xuYW1lcyhzY2FsZXMpIDwtIGMoIlBuIiwgIlFuIiwgIlNuIiwgIk1BRCIsICJTRCIsICJJUVIiLCAiWCIpDQoNCmRmIDwtIGRhdGEuZnJhbWUoc2NhbGVzKSAlPiUNCiAgZ2F0aGVyKGtleSwgdmFsdWUsIC1YKQ0KcCA8LSBnZ3Bsb3QoZGYsIGFlcyhYLCB2YWx1ZSwgY29sb3IgPSBrZXkpKSArIGdlb21fbGluZSgpDQpwbG90bHk6OmdncGxvdGx5KHApDQpgYGANCg0KDQoNCmBgYHtyfQ0KTiA8LSAxMDAwMA0KbXUgPC0gMg0Kc2lnbWEgPC0gMS8yDQpzY2FsZXMgPC0gbWF0cml4KG5jb2w9NikNCnBuIDwtIGMoKQ0KI1ggPSBybm9ybShOLCBtZWFuPW11LCBzZCA9IHNpZ21hKQ0KWF93ID0gcndlaWJ1bGwoTiwgc2hhcGU9bXUsIHNjYWxlPXNpZ21hKQ0KZm9yIChpIGluIHNlcSgwLjAxLCAwLjIsIGJ5PTAuMDEpKXsNCiAgWSA9IHJleHAoTippLCByYXRlPW11KQ0KICBYX1kgPSBjKFhfdywgWSkNCiAgcG4gPC0gYyhwbiwgUG4oWF9ZKSkNCiAgc2NhbGVzIDwtIHJiaW5kKHNjYWxlcywgYyhQbihYX1kpLCBRbihYX1kpLCBTbihYX1kpLCBtYWQoWF9ZKSwgc2QoWF9ZKSwgSVFSKFhfWSkpKQ0KfQ0Kc2NhbGVzIDwtIGNiaW5kKHNjYWxlcywgc2VxKDEsIDIxKSkNCmNvbG5hbWVzKHNjYWxlcykgPC0gYygiUG4iLCAiUW4iLCAiU24iLCAiTUFEIiwgIlNEIiwgIklRUiIsICJYIikNCg0KZGYgPC0gZGF0YS5mcmFtZShzY2FsZXMpICU+JQ0KICBnYXRoZXIoa2V5LCB2YWx1ZSwgLVgpDQpwIDwtIGdncGxvdChkZiwgYWVzKFgsIHZhbHVlLCBjb2xvciA9IGtleSkpICsgZ2VvbV9saW5lKCkNCnBsb3RseTo6Z2dwbG90bHkocCkNCmBgYA0KDQpgYGB7cn0NCk4gPC0gMTAwMDANCm11IDwtIDINCnNpZ21hIDwtIDEvMg0Kc2NhbGVzIDwtIG1hdHJpeChuY29sPTYpDQpwbiA8LSBjKCkNCiNYID0gcm5vcm0oTiwgbWVhbj1tdSwgc2QgPSBzaWdtYSkNClhfdyA9IHJ3ZWlidWxsKE4sIHNoYXBlPW11LCBzY2FsZT1zaWdtYSkNCmZvciAoaSBpbiBzZXEoMC4wMSwgMC4yLCBieT0wLjAxKSl7DQogIFkgPSByY2F1Y2h5KE4qaSwgbG9jYXRpb249bXUsIHNjYWxlPXNpZ21hLzEwMDApDQogIFhfWSA9IGMoWF93LCBZKQ0KICBwbiA8LSBjKHBuLCBQbihYX1kpKQ0KICBzY2FsZXMgPC0gcmJpbmQoc2NhbGVzLCBjKFBuKFhfWSksIFFuKFhfWSksIFNuKFhfWSksIG1hZChYX1kpLCBzZChYX1kpLCBJUVIoWF9ZKSkpDQp9DQpzY2FsZXMgPC0gY2JpbmQoc2NhbGVzLCBzZXEoMSwgMjEpKQ0KY29sbmFtZXMoc2NhbGVzKSA8LSBjKCJQbiIsICJRbiIsICJTbiIsICJNQUQiLCAiU0QiLCAiSVFSIiwgIlgiKQ0KDQpkZiA8LSBkYXRhLmZyYW1lKHNjYWxlcykgJT4lDQogIGdhdGhlcihrZXksIHZhbHVlLCAtWCkNCnAgPC0gZ2dwbG90KGRmLCBhZXMoWCwgdmFsdWUsIGNvbG9yID0ga2V5KSkgKyBnZW9tX2xpbmUoKQ0KcGxvdGx5OjpnZ3Bsb3RseShwKQ0KYGBgDQoNCg==